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1.0  Introduction 


This  is  the  report  for  fiscal  year  >992  on  the  project  “Aircraft  Wheel  Life  Assessment” 
funded  by  the  Landing  Gear  Systems  Section  of  the  Flight  Dynamics  Directorate  at  Wright- 
Patterson  Air  Force  Base. 

It  has  been  determined  that  an  important  part  of  the  wheel  life  assessment  problem  is  the  ac¬ 
curate  determination  of  the  tire-wheel  interface  pressure  distribution  under  various  loading  condi¬ 
tions.  A  combined  analytical/experimental  methodology  for  obtaining  this  pressure  distribution 
was  developed  at  the  University  of  Notre  Dame.  The  principle  analytical  tool  in  this  methodology 
is  the  finite  element  program  ANTWIL  (ANalysis  of  Tire-Wheel  Interface  Loads)  [Kandarpa  et 
al.,  1991)  which  recovers  the  pressure  distribution  given  a  number  of  experimental  strain  mea¬ 
surements  on  the  wheel.  The  major  activity  on  this  project  in  FY92  consisted  of  a  study  of  the  F- 
16  Block  30  and  the  Block  40  main  landing  gear  wheels  to  determine  the  optimal  number  and  lo¬ 
cation  of  the  strain  gages  for  subsequent  experiments.  Experiments  to  be  conducted  at  WPAFB 
will  record  strains  at  the  locations  specified,  and  these  data  will  be  used  to  determine  the  tire- 
wheel  interface  pressure  distributions  for  each  test  case.  Modifications  have  been  made  to 
ANTWIL  to  expand  its  capability.  The  ability  to  recover  a  circumferential  shear  traction  has  been 
incorporated  into  the  code.  Ongoing  work  and  planned  tasks  for  FY93  are  delineated. 

Landing  gear  wheel  failures  have  been  responsible  for  a  significant  amount  of  aircraft  dam¬ 
age.  Many  such  wheel  failures  have  been  catastrophic,  resulting  in  a  sudden  loss  of  tire  inflation 
pressure  and,  consequently,  severe  damage  to  primary  aircraft  structures.  An  FAA  study,  by  Du¬ 
mp  and  Brussat  [1985],  showed  that  more  than  30  percent  of  the  over  5,000  operational  failures 
reported  from  1970-1975  were  due  to  wheel  problems.  This  does  not  include  wheels  which  were 
removed  from  service  ahead  of  their  expected  service  life  because  of  cracks  found  during  routine 
tire  changes. 

The  Landing  Gear  Systems  Section,  WL/FIVMA,  at  Wright-Patterson  Air  Force  Base  initi¬ 
ated  an  inhouse  program  in  late  1986  to  investigate  experimental  and  analytical  methods  for 
wheel  life  estimation  and  verification  [Treanor  and  Carter,  1987;  Treanor  and  Boike,  1988].  As 
part  of  this  effort,  preliminary  studies  were  conducted  at  the  University  of  Notre  Dame  to  deter¬ 
mine  the  propagation  behavior  of  cracks  in  the  bead  seat  region  [Enneking,  1987;  Lawler,  et  al, 
1989].  While  the  results  of  the  effort  were  successful,  the  research  identified  a  need  to  determine 
the  load  distribution  at  the  tire-wheel  interface  and  the  distribution  of  time  to  crack  initiation  in 
order  to  accurately  predict  the  crack  propagation  life  of  a  wheel  in  the  bead  seat  region.  Subse¬ 
quent  wheels  mounted  with  radial  tires  have  failed  due  to  cracks  forming  in  other  parts  of  the 
wheel  {e.g.,  near  the  lightening  holes).  Therefore,  accurate  knowledge  of  the  stresses  at  arbitrary 
locations  in  the  wheel  is  necessary.  The  tire-wheel  interface  loads  are  also  necessary  for  such 
stress  determination. 
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The  University  of  Notre  Dame  has  completed  development  and  implementation  of  an  ex¬ 
perimental/analytical  methodology  which  will  allow  for  determination  of  the  loads  at  the  tire- 
wheel  interface.  A  finite  element  code,  ANTWIL,  has  been  developed  which  will  run  on  a  com¬ 
puter  workstation  in  the  WL/FIVMA  laboratory  and  can  directly  yield  the  tire-wheel  interface 
pressure  distribution  given  appropriate  strain  gage  measurements  [Kandarpa  et  al.,  1991].  The 
ANTWIL  program  employs  an  axisymmetric  finite  element  model  which  is  subjected  to  asym¬ 
metric  loading.  The  loading  is  represented  as  a  double  Fourier  series,  and  the  components  are  de¬ 
termined  by  a  least  squares  fit  using  the  experimentally  determined  strains. 

Below  we  review  the  work  completed  during  this  fiscal  year.  Following  this,  the  tasks  are 
discussed  that  have  been  identified  for  the  1993  fiscal  year  to  support  efforts  in  determining  air¬ 
craft  wheel  life. 
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2.0  Tire- Wheel  Interface  Load  Recovery 


2.1  Introduction 

One  of  the  research  efforts  in  FY  92  involved  finding  the  optimal  number  and  location  of 
strain  gages  necessary  to  capture  the  essential  features  of  the  tire-wheel  interface  pressure  distri¬ 
bution.  These  locations  were  subsequently  given  to  WL/FFVMA  laboratory  for  testing.  The  infor¬ 
mation  obtained  from  these  experiments  will  be  used  as  input  data  for  ANTWIL,  and  the  tire- 
wheel  interface  pressure  will  be  back-calculated  to  verify  the  proposed  methodology.  In  the  next 
few  paragraphs,  we  will  discuss  the  analytical  experiments  that  were  performed  for  determining 
the  optimal  number  and  location  of  strain  gages.  First,  the  load  recovery  algorithm  is  reviewed  for 
completeness. 

2.2  Theoretical  Development  of  Load  Recovery  Algorithm 

Consider  the  aircraft  wheel  depicted  in  Fig.  1.  The  geometry  of  the  wheel  and  strain  mea¬ 
surements  at  discrete  points  on  the  wheel  surface  are  assumed  to  be  available.  From  knowledge  of 
the  geometry  and  the  boundary  conditions  of  the  wheel,  an  axially  symmetric  finite  element  mod¬ 
el  is  developed  which  is  subjected  to  asymmetric  loading  (see  Appendix  A).  This  model  will  be 
used  to  back-calculate  the  load  distribution  at  the  tire-wheel  interface  (also  known  as  the  bead 
seat  region)  from  the  strain  measurements. 

Although  the  geometry  of  the  wheel,  in  general,  is  not  axially  symmetric  due  to  lightening 
holes,  brake  lugs,  etc.,  the  wheel  can  be  reasonably  approximated  as  an  axially  symmetric  body 
because  the  input  analysis  data  is  obtained  in  regions  removed  from  geometric  irregularities.  By 
using  this  geometric  simplification,  analysis  can  be  carried  out  on  a  two-dimensional  axisynunet- 
ric  domain  (r,  z)  rather  than  the  original  three-din^nsional  body,  thus,  greatly  reducing  the  com¬ 
putational  effort. 

In  the  proposed  problem  where  tire-wheel  interface  loads  are  to  be  computed,  the  loading  is 
not  only  asymmetric  but  also  varies  along  the  surface  (s)  of  the  body.  Herein,  we  use  a  cylindri¬ 
cal  coordinate  system  r,6^  (see  Fig.  2).  The  surface  load  is  expressed  as  a  function  of  0  and  5 
where  s,  the  surface  coordinate,  is  a  function  of  r  and  z.  The  loading  can  be  represented  as  a  Fou¬ 
rier  series  in  d  and  j.  In  the  algorithm  development  th^t  follows,  the  tractions  are  assumed  to  be 
normal  to  the  wheel  surface.  However,  in  the  theory  presented,  the  load  p  {6,  s)  can  represent  any 
loading  component,  including  surface  tractions.  The  load  p{d,s)  is,  thus,  given  by: 

pid,s)  =5^  [p^„  ( J)  cos  in  ff)  +  (s)  sin  (nd)  ] 
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nCURE  1:  Cross-Section  of  F-16  Block  30  Wheel  (showing  sectioned  area  only). 
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Further,  we  may  expand  the  coefficients  Pc„  (j) .  Ps„  (^)  in  Fourier  series  over  the  interval  of 
loading  defined  to  be  0  <  s  <  L.  Without  loss  of  generality  we  use  a  Fourier  cosine  series.  The 
final  series  representation  is  then  given  by: 


p{d,s)  =  2^2^{p^„^cos(nd)cos(m— ) +p,„^sin(n0)cos(m  — )  }.  (1) 

m  n 

Consider  a  single  loading  component  of  unit  amplitude: 

cos  (ntf)  cos  (m^) .  (2) 

Let  the  strain  response  due  to  the  above  load  be  (r,d,z)  ■  All  strains  calculated  in  this  paper 
are  based  on  a  linear  elastic  and  isotropic  material  response.  The  strain  tensor  at  a  point  due  to  this 
loading  can  be  decomposed  as  follows: 

z) cos (/i0) +€^„^(r, z)  sin (n^),  (3) 

where  and  are  the  cosine  and  sine  amplitudes.  If  the  load  amplitude  were  rather 
than  unity,  then  the  strain  response  would  be  Pc„„^cnm-  Similarly,  the  response  due  to  a  sine 
component  loading: 


sin  {n6)  cos  (m— ) 
in  the  ^-direction  would  be  written  as: 

=  «5nm  z)  COS  (nff)  +  (r,  z)  sin  (nO) . 


(4) 


(5) 


FIGURE  2.  A  figure  showing  the  r,  z  and  s  coordinates,  and  interface  pressure. 
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TTie  tensor  coefficients,  and  ^  ,  in  equations  (2)  and  (3)  have  the  structure. 


<  0  5^/2 

=  0^0. 


1^/2  0  < 


whereas  and  have  the  structure, 

cnm  snm  ’ 


I  0  2^*  0 

.0  y„  0 

The  total  response  due  to  the  loading  as  given  in  equation  (1),  by  superposition,  would  then  be: 


C  {rM  =  YL  ^Pcnm^cnm-^Psnm*-snm  > 

n  m 

=  XX  ^Pcnm  I'cnmCOS  {nS)  +€*„„sin  {nO)  ] 

n  m 

+  Psn«t«snmCOS(«^)  +e|„„sin(n^)]  } 

Since  strain  measurements  will  be  taken  on  the  surface  of  the  body,  an  expression  is  needed 
for  the  direct  strain  on  the  surface  in  a  specified  direction.  Referring  to  Fig.  3,  it  is  required  to  find 
the  strain  on  the  surface  of  the  body  of  revolution  at  a  point  with  surface  coordinate  s  and  in  a  di¬ 
rection  characterized  by  the  unit  vector  g,  defined  with  reference  to  Fig.  3  by  (see  Appendix  B): 

g  =  -/Lisin(^)e^+cos(<p)e^-»- Asin(^)e^  (8) 

where  and  A  are  direction  cosines  of  the  unit  tangential  vector  e,  and  <t>  is  the  angle  between 
the  strain  gage  placed  on  the  surface  of  the  body  and  the  circumferential  direction  e^. 

If  e  ( 0,j)  is  the  strain  tensor  of  equation  (8)  evaluated  at  the  point  of  interest  on  the  surface 
then  the  required  direct  strain  is  given  by: 

=  g^«(Mg  (9) 
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(a)  (b) 

FIGURE  3.  (a)  Axisynunetric  Body  (b)  IMrection  of  Vector  g  at  which  Strains  are  Computed. 

Thus,  we  can  write  the  series  expansion  for  this  direct  strain  as: 

€ie,s)  =  X2>  cnm  K««cos(n0)  +<„„sin(n^)l  + 

n  m  (10) 

P  snm  K««cos  (ntf)  +  4„sin  (n0)  ] 

where  the  subscript  g  is  dropped  for  simplicity,  but  it  should  be  kept  in  mind  that  all  the  strains  in 
equation  (10)  are  obtained  by  application  of  equation  (9)  to  the  component  strain  tensors  of  equa¬ 
tions  (6)  and  (7).  Also,  note  that  equation  (10)  may  be  reordered  as: 

=  XX  (^Pcnm<nm+Psnm^snm) 

n  m  (11) 

+P5««4«)  Sin  {nd) 

The  direct  strains  are  obtained  experimentally  at  J  locations  and  are  denoted  e*  (sj) ; 
j  =  1, .. .,  7  (in  the  direction  defined  by  g  as  above).  They  may  be  decomposed  into  Fourier  co¬ 
sine  and  sine  components  and  written  as: 

iSj)  =  X  ^  ^*sn  sin  (n  0)  }  (12) 

n 

Now  comparing  the  cosine  and  sine  terms  of  equations  (11)  and  (12)  yields. 
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(13) 


M 

m  =  0 

and: 


M 

^Pcnm^cnm^^j^  Psnm^snm^^j^  ^  ~ 

m  -  0 

where  the  series  has  been  truncated  after  M  terms. 

The  above  equations  can  be  concisely  written  in  matrix  notation  as: 


where: 


E-,P„  =  €* 


nj' 
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(14) 


(15) 


4i^)  ••• 

(16) 

4i  4i  (^>)  - 

••• 

1*«  “  \Pcn\  Psn\  *• 

PcnM 

(17) 

These  equations  are  assembled  together  for  all  J  experimental  points  resulting  in: 


(18) 


E„P„  =  €*  (19) 

where: 

K  =  [E.i  -  sj' 

<  =  [<,  -‘Sj]'  (20) 

Finally,  premultiplying  equation  (19)  by  results  in  the  least  squares  equation, 

(EX)P,  =  (21) 


s 


Solution  of  equation  (21)  yields  the  amplitudes  of  the  Fourier  components,  p^.  This  process  is 
repeated  for  all  n  =  0, iV  in  the  ^-direction.  The  amplitudes  are  substituted  into  equation 
(1)  and  the  load  distribution  on  the  contact  surface  is  produced. 

23  Numerical  Experiments 

Sample  numerical  experiments  were  conducted  in  a  three  step  sequence.  First,  a  known  an¬ 
alytical  load  (see  Fig.  4)  is  applied  on  the  surf^e  of  the  wheel  (around  the  bead  seat  region).  By 
performing  a  finite  element  analysis,  strain  readings  are  obtained  at  J  locations  on  the  surface. 
Usually  the  strains  were  computed  at  the  same  locations  as  the  experimental  strains  gages  and  for 
this  reason  are  called  the  experimental  strains.  These  strains  are  decomposed  in  their  sine  and  co¬ 
sine  components  and  stored  in  a  vector  of  length  27,  e*  ,  where  n  corresponds  to  the  harmonic 
number  in  the  circumferential  direction.  Note  that  in  an  actual  analysis,  the  first  step  is  done  ex¬ 
perimentally.  Next,  the  same  analytical  model  is  subjected  to  M  unit  harmonic  loads  in  the  j -di¬ 
rection,  for  a  specific  n,  and  strains  are  determined  at  the  same  locations  as  the  experimental 
strains.  These  analytical  strains  are  then  stored  in  a  matrix  of  dimension  2JxM  where  each 
column  of  length  27  corresponds  to  the  strains  due  to  the  m‘**  Fourier  load.  Finally,  the  Fourier 
amplitudes  p„  are  obtained  by  solution  of  the  least  squares  equation  (20).  These  amplitudes  are 
used  for  calculating  the  interface  pressure  which  is  then  compared  with  the  applied  load  for  con¬ 
vergence  and  accuracy.  Further  details  of  the  above  procedure  can  be  found  in  Kandarpa  et  al. 
(1991). 

Using  the  above  procedure,  sample  numerical  experiments  were  conducted  by  using  differ¬ 
ent  numbers  of  strain  gages.  These  results  were  further  compared  among  themselves  and  the  best 
configuration  was  chosen.  Details  of  the  analytical  experiments  are  summarized  below. 


FIGURE  4.  Hypothetical  Load  on  Bead  Seat  Regkm 
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4b  diff  of  areas—  1 .62 1 
L.2  norm—  13.9 


FIGURE  5.  Comparison  of  Anaiytkal  and  Recovered  Load.  (»• 
recovered  load,  —  applied  load;  Rosettes:  4  out,  4  in,  2  flange) 


23.1  Block  30  Wheel 

From  research  performed  earlier,  it  was  determined  that  the  best  results  were  achieved  when 
there  was  a  significant  number  of  strain  gages  in  the  loading  area.  It  was  found  that  eight  strain 
gage  rosettes  (total,  inboard  and  outboard  sides  of  the  wheel)  in  the  loading  area  provided  ade¬ 
quate  results.  Two  additional  rosettes  were  placed  at  tlM  base  of  the  flange  to  obtain  a  better  recov¬ 
ery  of  the  load.  The  need  for  the  gages  on  the  base  of  the  flange  was  also  determined  from 
previous  experiments. 

The  number  of  strain  gage  rosettes  in  the  bead  seat  region  were  evenly  placed  along  the 
flange  on  both  the  inside  and  outside  surfaces.  Four  rosettes  were  placed  on  the  inside  surface  and 
four  were  placed  on  the  outside,  along  with  the  two  gages  on  the  base  of  the  flange,  for  a  total  of 
10  gages.  Figure  5  shows  the  analytically  ^plied  load  as  a  dashed  line  and  the  recovered  load  as 
a  solid  line.  This  configuration  was  used  as  a  base  of  coiiq;>arison  as  the  interior  rosettes  were  re¬ 
moved  one  at  a  time. 

Three  measures  of  error  are  given  on  each  graph.  The  L^-norm  (denoted  the  infinity  norm) 
is  the  absolute  value  of  the  maximum  difference  between  the  applied  load  (s)  and  the  recov¬ 
ered  load  p*  ,  Le., 


L^-norm  =  max  |  p^{s)  -p*{s)  \  .  (22) 

J 

The  L2-norm  is  the  square  root  of  the  square  of  the  differences  integrated  over  the  loading  region, 
Le., 
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L2-norm  =  {J|  p,(j)  -p^is) 


(23) 


The  third  error  measure,  closely  related  to  the  L,  -  norm,  is  called  here  the  percent  diff  of  areas, 
and  is  defined  as: 

fl  P«(^)  I 

%  diff  of  areas  =  = - p - x  IfX).  (24) 

]p^{5)  ds 

The  results  depicted  in  Fig.  6  were  obtained  from  the  equal  spacing  of  three  rosettes  on  the 
inside  curve  for  a  total  of  nine  rosettes.  All  error  measures  show  this  configuration  to  be  marginal¬ 
ly  better  than  the  previous  configuration. 

The  number  of  equally  spaced  rosettes  on  the  inside  of  the  load  area  was  then  reduced  to 
two.  Again,  the  desired  result  was  to  decrease  the  number  of  gages  on  the  inside  surface  without 
losing  accuracy.  The  results  are  displayed  in  Fig.  7.  The  errors  are  close  to  those  of  the  original 
configuration  except  for  the  L2-norm  which  is  higher. 

As  shown  in  the  previous  figures,  the  optimal  configuration  in  the  loading  area  includes  four 
rosettes  on  the  outside  of  the  load  area  and  three  rosettes  on  the  inside  of  the  load  area.  The  next 
step  was  to  determine  the  variations  in  the  results  that  could  be  caused  by  adding  a  rosette  to  the 
area  at  the  base  of  the  flange.  The  configuration  also  included  four  rosettes  on  the  outside  of  the 


FIGURE  6.  CompariwMi  of  Analytical  and  Recovered  Load.  (— 
recovered  load,  —  applied  load;  Rosettes:  4  out,  3  in,  2  flai^) 
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FIGURE  7.  C^omparkoii  of  Analytkal  and  Recovered  Load.  (-> 
recovered  kMd,  —  appUed  load;  Roeettes:  4  ont,  2  fait  2  fange) 


load  area,  three  on  the  inside  of  the  load  area,  and  the  two  rosettes  at  the  base  of  the  flange.  Figure 
8  illustrates  the  results  for  this  configuration. 

Comparing  these  results  with  the  results  from  the  configuration  without  the  extra  rosette  on 
the  base  of  the  flange,  it  can  be  seen  that  the  extra  rosette  provides  redundant  information.  The  re¬ 
sulting  recovered  load  has  a  slightly  higher  error  for  two  of  the  three  error  norms  with  the  extra 
rosette  on  the  base  of  the  flange. 

The  next  run  was  to  determine  if  adding  a  rosette  on  the  base  of  the  flange  would  improve 
the  results  from  the  configuration  with  only  two  rosettes  on  the  inside  and  four  on  the  outside  of 
the  loading  surface.  By  adding  the  gages,  the  results  depicted  in  Fig.  9  were  recovered.  The  results 
of  this  configuration,  relative  to  the  result  of  the  corresponding  configuration  without  the  extra 
gage.  Fig.  7,  are  approximately  the  same.  The  third  rosette  in  the  middle  of  the  flange  i»rvides  lit¬ 
tle  additional  information  and  is  not  needed. 

The  next  experiments  were  done  to  determine  the  best  placement  of  the  two  rosettes  on  die 
base  of  the  flange.  In  the  previous  trials,  the  base  rosettes  were  in  a  set  place  and  were  not  moved. 
In  Fig.  10,  the  base  rosettes  were  moved  one  finite  element  towards  the  left.  The  load  recovered  is 
essentially  equivalent  to  the  load  founu  in  Fig.  6.  However,  moving  the  rosettes  another  element 
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Load(psi)  Load(psi) 


FIGURE  8.  Comparison  of  Analytkal  and  Recovered  Lioad.  (— 
recovered  load,  —  aisled  load;  Rosettes:  4  out,  3  In,  3  flange) 


FIGURE  9>  Comjparison  of  Analytical  and  Recovered  Load.  (—  recovered 
load,  —  applied  load;  Rosettes:  4  oat,  2  in,  3  flange) 
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to  the  left  produced  significantly  poorer  results.  It,  therefore,  was  concluded  that  the  configuration 
in  Fig.  6  was  the  best  configuration  of  rosettes  to  adequately  recover  the  load. 

The  previously  presented  sample  problems  have  been  analyzed  using  the  zeroth  harmonic 
in  the  circumferential  direction.  The  results  of  other  harmonics  will  now  be  presented  and  dis¬ 
cussed. 

The  results  for  the  first,  second,  and  third  harmonics  are  in  Figs.  11,  12,  and  13,  respective¬ 
ly.  The  results  for  all  three  harmonics  are  comparable  to  the  zeroth  harmonic  results  for  the  same 
configuration  of  strain  gages. 

Finally,  to  determine  the  robustness  of  the  configuration  of  the  strain  gages,  a  different  type 
of  load  was  applied.  In  this  case,  the  applied  load  was  the  mirror  image  of  the  previous  one.  This 
load  is  shown  as  the  solid  line  in  Fig.  14.  The  same  configuration  of  strain  gages  as  used  in  Fig.  6 
is  used  for  the  analysis.  The  analysis  was  performed  for  the  zeroth  mode,  and  the  result  is  repre¬ 
sented  by  the  dashed  line  in  Fig.  14.  Although  tltt  results  are  not  quite  as  accurate  as  those  in  Fig. 
10,  it  can  be  concluded  that  even  with  a  severe  modification  of  the  loading,  this  strain  gage  config¬ 
uration  allows  an  accurate  recovery  of  the  load. 


FIGURE  10.  Comparison  of  Analytical  and  Recovered  Load  Varying  Location  on 
Flange.  ( —  recovered  load,  —  applied  load;  Rosettes:  4  ont,  3  in,  2  flange) 
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Lo>d(psi)  Loid(iKi) 


FIGURE  11.  Results  of  ANTWIL  for  Pint  Mode  for  Loadf  Applied  and  Recovered  io  the 
s(r4)  Direction.  (-~  recovered  loud,  —  applied  load;  Rosettes:  4  out,  3  in,  2  lange) 


FIGURE  12.  Residts  of  ANTWIL  for  Second  Mode  for  Loads  Apidkd  and  Recovered  in  tbe 
s(r,z)  Direction.  ( — recovered  load,  —  applied  load;  Rosettes:  4  out,  3  in,  2  flange) 
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2.3.2  Block  30  Wheel:  Optimal  Locations 

The  various  sample  problems  that  were  carried  out  on  the  Block  30  wheel  were  helpful  in 
locating  the  key  positions  for  the  strain  gages  in  the  bead  seat  region  of  the  wheel.  These  positions 
provide  the  optimal  strain  information  allowing  the  algorithm  to  best  capture  the  tire-wheel  inter¬ 
face  pressure  distribution.  Figure  IS  and  Table  1  illustrate  these  strain  gage  locations. 


Table  1:  Gage  Positions*  (F-16  Block  30  Wheel) 


Gage(s) 

Z  -Position 
(inches) 

Side  of 
Flange 

Orientation** 

Distance 

Along 

Flange 

1 

0.6 

Tire*" 

Rosette 

2.38  inches 

2 

0.98 

Tire 

Rosette 

3 

1.6 

Tire 

Rosette 

4 

1.68 

Outside** 

Rosette 

5 

1.11 

Outside 

Rosette 

6 

0.75 

Outside 

Rosette 

7 

0.39 

Outside 

Rosette 

8 

3.35 

Outside 

Circumferential 

9 

3.35 

Tire 

Circumferential 

a.  The  gages  are  to  be  placed  in  line  with  a  large  lightening  hole  (away  from  brake  lug). 

b.  If  the  orientation  is  a  rosette,  place  a  rosette  with  middle  leg  in  the  z  direction  from  in¬ 
board  to  outboard,  i.e.,  from  left  to  right  of  Fig.  IS. 

c.  This  side  of  the  flange  refers  to  the  side  where  the  bead  seat  and  the  tire  are  placed. 

d.  This  side  of  the  flange  refers  to  the  side  opposite  from  the  tire. 


2.3.3  Block  40  Wheel 

The  blue-print  of  the  Block  40  wheel  supplied  by  the  Landing  Gear  Systems  Section  at 
Wright-Patterson  Air  Force  Base  was  used  in  generating  an  axisymmetric,  eight-node  finite  ele¬ 
ment  mesh  using  Patran  (see  Fig.  16).  This  mesh  was  used  for  the  analysis  of  load  recovery  at  the 
tire-wheel  interface.  The  Block  30  results  were  then  used  as  a  guide  for  the  initial  strain  gage 
placement  on  the  Block  40  wheel.  Sample  problems  for  the  Block  40  were  analyzed  by  using  the 
zeroth,  first,  second,  and  third  harmonics  in  the  circumferential  direction.  The  mirror  image  of  the 
current  load  was  also  applied  and  recovered.  These  results  are  presented  in  this  report. 
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FIGURE  15.  Locations  of  Strain  Gafcs  for  F-16  Block  30  Wbcd. 


Using  the  Block  30  wheel  as  a  model,  nine  strain  gage  rosettes  were  used  in  the  Block  40 
test.  Seven  gages  were  placed  in  the  bead  seat  region,  three  on  the  tire  side,  four  on  the  outside. 
They  were  spaced  to  achieve  an  even  distribution  along  the  loading  area.  The  remaining  two  test 
points  were  at  the  base  of  the  flange,  one  on  the  tire  side,  one  on  the  outside.  This  configuration 
returned  the  load  very  accurately.  Figure  17  shows  the  applied  load  and  the  recovered  load  along 
with  the  error  measures  employed  to  assess  the  effectiveness  of  the  recovery.  The  two  flange  area 
gages  were  moved  toward  the  bead  seat  area  to  determine  if  the  results  could  be  improved,  but 
this  was  found  to  have  a  detrimental  effect. 

To  see  if  the  two  gages  in  tlK  flange  area  were  needed  to  obtain  satisfactory  results  on  the 
Block  40  wheel,  the  gages  in  the  flange  area  were  removed  entirely,  leaving  only  the  seven  in  the 
bead  seat  area.  This  change  also  had  a  deleterious  effect.  To  make  up  for  the  information  lost 
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when  the  two  flange  gages  were  removed,  an  eighth  rosette  was  added  on  the  outside  of  the  bead 
seat  area  and  the  eight  gages  were  again  positioned  to  obtain  an  even  spatial  distribution.  These 
results  were  not  quite  as  good  as  the  first  set  with  the  flange  gages  (see  Fig.  18).  However,  this 
configuration  of  eight  gages  was  superior  to  the  nine  gage  configuration  for  all  the  higher  harmon¬ 
ics  tested.  The  recovery  for  harmonics  one,  two,  and  three  are  shown  in  Figs.  19,  20,  and  21.  A 
ninth  gage  was  added  to  the  bead  seat  area  to  determine  if  more  information  could  be  gained,  but 
this  proved  ineffective. 

Thus,  best  configuration  is  obtained  using  eight  gages,  all  in  the  bead  seat  area.  This  config¬ 
uration  was  also  used  to  successfully  recover  a  mirror  image  of  the  modified  load  for  the  zeroth, 
first,  second,  and  third  modes.  See  Figs.  22  through  25. 

Based  on  the  results  of  these  studies.  Fig.  16  and  Table  2  indicate  the  recommended  place¬ 
ment  of  the  eight  gages  required  for  satisfactory  tire-wheel  interface  pressure  recovery. 


Table  2:  Gage  Positions  (F-16  Block  40  Wheel) 


Gage(s) 

Z  -Position 
(inches) 

Side  of 
Flange 

Orientation® 

Distance 

Along 

Flange 

1 

0.73 

Rosette 

0.94  inches 

2 

0.48 

Outside‘s 

Rosette 

3 

0.87 

Tire 

Rosette 

4 

0.72 

Outside 

Rosette 

5 

1.06 

Tire 

Rosette 

6 

0.98 

Outside 

Rosette 

7 

1.61 

Tire 

Rosette 

8 

1.59 

Outside 

Rosette 

9 

2.89 

Tire 

Circumferential 

10 

2.89 

Outside 

Circumferential 

a.  If  the  orientation  is  rosette,  place  a  rosette  at  this  position  with  the  middle  leg  in  the  z 
direction  from  inboard  to  outboard,  i.e.,  from  left  to  ri^t  of  Fig  16. 

b.  This  side  of  the  flange  refers  to  the  side  where  the  bead  seat  and  tire  are  placed. 

c.  This  side  of  the  flange  refers  to  the  side  opposite  from  the  tire. 
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FIGURE  16.  Locations  of  the  Optimal  Number  of  Strain  Gages  for  the  F>16  Block  40  Wheel. 


2A  Shear  TVactions 

The  analytical/experimental  algorithm  has  been  extended  to  include  determination  of  tire- 
wheel  interface  shearing  tractions,  and  the  ANTWIL  finite  element  code  was  modified  to  include 
this  effect.  Complex  wheel  loadings  which  occur  in  braking,  spin-up,  cornering,  etc.,  produce 
both  normal  and  shear  tractions  at  the  tire-wheel  interface.  The  effect  of  the  shear  tractions  is 
thought  to  be  important  to  durability  and  damage  tolerance  of  the  wheel. 

The  previous  coding  in  ANTWIL  was  based  on  calculating  nodal  loads  by  using  normal 
components  of  applied  surface  tractions.  The  code  has  been  modified  to  allow  surface  tractions  of 
arbitrary  orientation.  This  permits  surface  tractions  with  arbitrary  normal  and  shear  components. 
It  should  be  noted  that  the  algorithm  given  in  Section  2.2  applies  to  the  shear  traction  problem  if 
the  pressure  p{d,s)  is  taken  to  be  the  shear  traction.  Of  course,  the  finite  element  load  vector 
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Load  (psi)  (P*') 


FIGURE  17.  Results  of  ANTWIL  for  Nine  Gage  Configuration  ( —  recovered  load, 
—  applied  load;  Rosettes:  4  out,  3  in,  2  flange) 


FIGURE  18.  Results  of  ANTWIL  for  Zeroth  Harmonic  for  Loads  Applied  and 
Recovered  in  the  s(r,z)  Direction:  Eight  Gage  Configuration  (•••  recoveiVd  load,  — 
applied  lom§;  Rosettes:  4  out,  4  in.) 
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FIGURE  19.  Results  of  ANTWIL  for  First  Harmoak  for  Loads  Applied  and  Recovered  in  the  s(r^) 
Direction:  Eight  Gage  Configuration  ( —  recovered  load,  —  applied  load;  Rosettes:  4  out,  4  iiL) 


FIGURE  20.  Results  of  ANTWIL  for  Second  Harmonic  for  Loads  Applied  and  Recovered  in  the 
i(r,z)  DirMtion:  Eight  Gage  Configuration  ( —  recovered  load,  —  applM  load;  Rosettes:  4  out,  4  in.) 
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FIGURE  22.  Results  of  ANTWIL  for  Zeroth  Harmonic  Using  a  Mirror  Image  Load: 
Eight  Gage  Configuration  ( —  recovered  load,  —  applied  load;  Rosettes:  4  out,  4  in.) 
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Load(psi) 


FIGURE  23.  Results  of  ANTWIL  for  First  Harmonic  Using  a  Mirror  Image  Load: 
Eight  Gage  Configuration  ( —  recovered  load,  —  applied  load;  Rosettes:  4  out,  4  in.) 


FIGURE  24.  Results  of  ANTWIL  for  the  Second  Harmonic  Using  a  Mirror  Image  Load: 
Eight  Gage  Configuration  ( —  recovered  load,  —  applied  load;  Rosettes:  4  out,  4  in.) 
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would  be  assembled  accordingly.  Several  sample  problems  of  shear  traction  loading  have  been 
analyzed  to  verify  the  program  modifications. 

2.4.1  Numerical  Experiments. 

This  research  was  done  to  determine  if  the  proposed  number  and  locations  of  strain  gages 
were  adequate  to  recover  tire-wheel  interface  shear  loads.  Trials  were  first  done  using  a  simple 
cylinder  in  which  the  exact  answer  was  known.  These  experiments  allowed  for  a  check  of  the 
ability  of  ANTWIL  to  recover  a  load  in  shear  and  have  the  model  react  in  an  expected  way.  The 
cylinder  was  modeled  using  a  fine  and  a  course  mesh,  both  of  which  produced  similar  results,  in¬ 
dicating  an  adequate  degree  of  numerical  convergence  had  been  attained.  Following  these  results, 
trials  were  run  on  the  Block  30  wheel. 

Experiments  conducted  on  the  wheel  consisted  of  various  loadings,  including  a  constant 
loading  and  an  assumed  variable  load.  ANTWU.  was  able  to  accurately  recover  these  loads.  The 
load  employed  had  the  same  profile  that  was  applied  normal  to  the  surface  in  previous  trials.  This 
load  is  depicted  in  Fig.  26.  The  recovered  load  is  superimposed  over  the  ^plied  load  and  error 
norms  are  given.The  figure  illustrates  the  degree  to  which  the  program  was  able  to  recover  the  ap¬ 
plied  load.  These  results  were  not  as  accurate  as  those  achieved  in  the  normal  loading.  However, 
the  shear  loads  are  usually  considered  of  secondary  importance  in  the  determination  of  the  wheel 
life. 
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FIGURE  26.  Rcsaitg  of  ANTWIL:  Recovery  of  Shear  T^actkMM 
(—  racoveicd  load,  —  applied  load;  Roeettes:  4  oat,  3  in,  2 100(0) 


3.0  Crack  Initiation:  Literature  Review 


3.1  Introduction 

It  is  well  known  that  the  combination  of  cyclic  stresses  due  to  wheel  rotation  and  corrosive 
environmental  influences  can  cause  surface  cracks  to  form  which  further  enhance  stress  concen¬ 
trations  leading  to  crack  growth  and  ultimate  failure  of  the  wheel.  While  current  durability  and 
damage  tolerance  approaches  neglect  the  initiation  life  of  fatigue  cracks,  this  assumption  appears 
to  be  too  conservative  for  wheel  life  predictions  and  needs  to  be  reevaluated  in  the  context  of 
overall  life  prediction  methodology. 

A  large  number  of  metallurgical  observations  and  numerous  experiments  have  been  carried 
out  to  understand  the  basis  of  crack  nucleation.  The  models  formed  to  understand  the  process  of 
crack  nucleation  are  based  on  a  wide  variety  of  theories  which  are  both  analytical  and  phenome¬ 
nological.  In  the  following  sections,  an  effort  has  been  made  to  separate  these  models  into  a  few 
important  groups  which  share  a  common  basis  or  approach.  This  review  is  based  primarily  on  a 
literature  survey  performed  at  WPAFB  and  with  the  assistance  of  WPAFB  personnel.  The  models 
surveyed  can  be  divided  into  three  principal  categories:  slip-band  theories,  continuum  damage 
mechanics  theories,  and  stochastic  models.  A  brief  review  of  this  literature  including  a  brief  histo¬ 
ry  of  the  evolution  of  the  models  is  presented. 

3.2  Models  Based  on  Slip  Band  Theory 

According  to  this  theory,  the  site  of  crack  initiation  is  often  associated  with  bands  of  intense 

slip  which  are  often  called  persistent  slip  bands  (PSBs).  These  bands  are  formed  due  to  the  pres¬ 
ence  of  dislocations  in  the  material  lattice.  The  dislocations  are  caused  either  due  to  exclusions  or 
intrusions  of  foreign  particles  in  the  structural  matrix. 

The  first  serious  study  of  crack  nucleation  based  on  slip  band  theory  was  conducted  by  Th¬ 
ompson,  et  al.  [1953].  The  experimental  studies  performed  on  copper  specimens  by  this  group 
showed  that  cracks  often  nucleated  at  the  surface  of  the  material  where  slip  bands  were  present.  In 
fact,  the  experiments  showed  enhanced  performance  for  the  specimens  in  which  these  persistent 
slip  bands  were  removed  physically,  by  electro-polishing,  during  the  course  of  the  experiments. 
There  has  followed  a  continuous  research  effort  to  investigate  the  surface  nucleation  of  these  slip 
bands. 

Some  of  the  earlier  models  were  based  on  Mott’s  assumption  [Mott,  1958],  that  the  disloca¬ 
tions  moved  in  different  paths  in  slip  bands  under  cyclic  loading.  But  these  models  failed  to  ex¬ 
plain  the  systematic  buildup  of  plastic  strain  during  a  cyclic  test.  Lin  and  Ito  [1969]  proposed  that 
one  part  of  the  slip  band  slides  in  the  forward  direction  and  the  other  part  slides  in  the  reverse  di¬ 
rection.  To  make  this  mode  of  deformation  possible,  they  assumed  a  special  internal  stress  field. 
This  theory  was  largely  ignored  because  it  was  considered  that,  in  reality,  initial  stresses  were  not 
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present.  Since  that  time,  however,  researchers  (e.g..  Brown  and  Ogin  [1985]  and  Laird  [1979]) 
have  established  that  internal  stresses  are  indeed  present.  Although  the  presence  of  internal  stress¬ 
es  in  the  PSB  has  been  established,  there  is  still  some  debate  as  to  whether  they  remain  constant  or 
change  during  cyclic  loading,  and  this  has  been  a  major  research  area  for  some  time.  Lin  [1985] 
used  the  theory  of  internal  stresses  in  slip  bands  to  estimate  the  damage  accumulated  during  cyclic 
loading. 

A  model  based  on  slip  bands  was  also  proposed  by  Tanaka  [1981].  This  model  differed 
from  previous  models  by  assuming  that  damage  was  caused  by  plastic  deformation  within  the  slip 
band  and  was  modeled  by  two  adjacent  layers  of  dislocation  pile  ups. 

Efforts  were  also  made  to  investigate  the  structure  and  properties  of  PSBs  and  their  in¬ 
volvement  in  crack  initiation.  Mughrabi  [1980],  and  later  Brown  [1977,  1980],  established  that 
PSBs  are  bands  of  concentrated  slips  which  essentially  carry  all  the  plastic  strain  in  low  strain  am¬ 
plitude  fatigue  test.  Some  studies  also  tried  to  address  the  issues  of  the  slowing  of  crack  nucle- 
ation  at  grain  boundaries  and  the  effects  of  grain  size  on  the  threshold  for  fatigue  crack  growth 
[Ravindran  and  Dwarkadasa,  1991].  Yoo  [1980]  conducted  a  study  to  investigate  twinning  and 
cleavage  fracture  using  a  dislocation  model.  The  model  was  based  on  the  anisotropic  elasticity 
theory  of  dislocations.  High  temperature  effects  in  PSBs  have  also  been  investigated.  lizuka  and 
Tanaka  [1991]  investigated  temperature  effects  on  stainless  steel  at  973  K  for  low  cycle  fatigue 
tests. 

A  different  approach  was  taken  by  Morris  et  al.  [1988].  The  nnodel  assumed  that  slip  band 
nucleation  and  the  strain  hardening  parameters  are  functions  of  local  plastic  flow.  The  functions 
incorporated  a  crack  growth  equation  that  was  modified  for  plastic  deformation  to  calculate  the 
cycles  required  for  a  microcrack  to  reach  the  first  subsurface  grain  boundary  when  propagating  in 
a  region  of  localized  plastic  flow. 

3.3  Continuum  Damage  Mechanics  Models 

Between  the  microscale  ( 10"^  - 10"^  mm)  and  structural  ( 10^  -  10^  mm)  scale,  there  exists  a 

macroscale  level.  The  continuum  damage  mechanics  approach  deals  within  that  macroscale  by 
defining  a  damage  variable  as  an  effective  surface  density  of  cracks  or  cavity  intersections  with  a 
plane.  Kachanov  [1958]  used  an  effective  mean  stress  as  a  damage  variable  to  model  creep  rup¬ 
ture.  This  work  was  the  starting  point  of  other  continuum  damage  models  extended  to  include  dis¬ 
sipation  and  low  cycle  fatigue  in  metals  by  Lemaitre  [1971]  and  coupling  between  damage  and 
creep  by  Leckie  [1974]. 

Some  of  the  earliest  work  on  a  continuum  damage  mechanics  model  for  fatigue  cracks  was 
carried  out  by  McClintock  [1968]  and  Rice  and  Tracey  [1969].  These  models  considered  defects, 
i.e.,  the  growth  of  cavities  and  microcracks  by  analyzing  their  geometry  in  a  continuous  matrix 
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using  principles  of  continuum  mechanics.  TlMse  models  explained  the  phenomenon  of  crack  initi¬ 
ation  by  means  of  a  micro-mechanics  analysis.  However,  it  was  recognized  that  microscale  diffi¬ 
culties  arise  when  the  analysis  has  to  be  applied  at  a  structural  scale  to  predict  ductile  failures. 

Chaboche  [1978]  provided  a  continuum  thermodynamic  foundation  to  justify  the  continu¬ 
um  damage  mechanics  theory.  He  also  studied  structural  behavior  by  introducing  internal  vari¬ 
ables  to  develop  a  number  of  phenomenological  laws.  The  formulation  adopted  three  types  of 
internal  variables,  namely,  an  isotropic  hardening  variable,  a  kinematic  hardening  variable  and  a 
damage  variable.  The  model  emphasized  the  relationship  between  damage  and  behavior  in  defor¬ 
mation.  Lemaitre  [1985]  formulated  a  model  that  used  effective  stress  and  the  laws  of  thermody¬ 
namics  to  assess  damage.  The  damage  in  this  case  was  linear  with  equivalent  strain  and  showed  a 
large  influence  of  triaxial  effects  through  the  equivalent  stress.  However,  this  model  was  limited 
by  the  hypothesis  of  isotropy  of  damage  and  isotropy  of  plasticity. 

A  continuum  mechanics  model  by  Tanaka  et  ai  [1988]  incorporated  a  recovery  effect  by 
diffusion  of  atoms  to  explain  the  initiation  of  grain-boundary  wedge-type  cracks  in  high-tempera- 
ture  fatigue.  They  assumed  that  the  microscopic  direction  of  deformation  is  periodically  changed 
in  fatigue  and  that  crack  initiation  was  dependent  on  the  deformation  history  of  the  grain  bound¬ 
ary. 


3.3.1  Stochastic  Models 

Laboratory  experiments  show  that  even  under  identical  laboratory  conditions  and  uniform 
cyclic  loading,  macroscopically  identical  metal  alloy  specimens  exhibit  large  statistical  variations 
in  the  number  of  fatigue  cycles  needed  to  create  a  macroscopic  crack.  It  can  be  argued,  therefore, 
that  fatigue  failure  theory  is  properly  cast  as  a  problem  in  probabilistic  structural  mechanics. 

Stochastic  models  in  fatigue  life  prediction  have  appeared  in  the  literature  for  some  tiiiM, 
primarily  devoted  to  the  prediction  of  crack  growth  (see  Sobczyk  and  Spencer  [1992]).  Cox  and 
Morris  [1987]  formulated  a  probabilistic  model  which  attempts  to  link  stochastic  microstructure 
and  the  statistics  of  measured  growth  rates.  The  model  is  formulated  as  a  semiMarkov  processes. 
The  underlying  Markov  process  describes  the  evolution  of  a  growth  control  variable  as  an  explicit 
function  of  crack  length.  Elapsed  fatigue  cycles  and  the  distribution  of  times  of  failure  are  calcu¬ 
lated  by  invoking  an  empirical  or  postulated  law  of  growth  rate.  This  law  is  either  a  deterministic 
or  probabilistic  relationship  between  the  growth  control  variable  and  the  crack  velocity.  The  mod¬ 
el  also  facilitated  the  use  of  spectral  loading. 

In  another  paper  (Cox,  et  al.  [1987])  assumed  that  the  stochastic  nature  of  the  short  surface 
-  fatigue  cracks  (on  the  order  of  the  grain  size)  under  uniform  loading  may  be  attributed  to  their 
interaction  with  the  stochastic  microstructure.  In  this  paper,  a  computational  scheme  was  pro¬ 
posed  to  describe  the  statistical  effects  of  the  temporary  arrest  of  short  cracks  by  grain  boundaries, 
and  the  modification  of  rates  of  propagation  across  grains  by  various  effects  of  grain  boundary 
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constraint.  However,  this  model  was  restricted  to  Mode  I  (tensile  opening)  transgranular  crack 
growth  under  constant  amplitude  cyclic  loading. 

The  model  by  Graef  [1981],  based  on  slip  dislocations,  used  probabilistic  theory  to  compute 
the  influence  of  the  microscopic  inhomogeneity  of  plastic  deformation  on  initiation  of  fatigue 
cracks.  A  connection  between  the  probability  of  the  occurrence  of  high  surf£a:e  steps  and  the  num¬ 
ber  of  stress  cycles  up  to  crack  initiation  was  noted.  He  showed  that  it  was  possible  to  calculate 
the  probability  of  occurrence  of  initial  cracks  and  to  correlate  the  results  of  investigations  on  the 
first  stage  of  crack  propagation  in  various  alloys. 

33^  Simunary 

The  models  presented  in  the  above  summary  review  several  different  theories  used  for  un¬ 
derstanding  crack  growth  and  crack  nucleation  in  metals  It  is  clear  from  the  literature  that  a  satis¬ 
factory  approach  to  the  problem  of  the  prediction  of  crack  initiation  is  still  forthcoming. 
Additionally,  the  fertile  subject  of  probabilistic  mechanics  seems  not  yet  to  have  yielded  suffi¬ 
ciently  useful  results  to  allow  its  direct  application  to  crack  initiation  in  aircraft  wheels.  Also  note 
that  the  models  presented  in  this  review  do  not  include  those  that  are  based  on  continuum  fracture 
mechanics,  that  is,  those  that  use  macroscopically  observable  parameters  to  predict  crack  initia¬ 
tion.  These  models  were  largely  derived  for  crack  propagation  and  were  later  extrapolated  to  initi¬ 
ation  problems. 
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4.0  Residual  Stresses:  Literature  Review 


4.1  Introduction 

Mechanical  components  are  often  subjected  to  surface  processing  treatments  to  improve 
mechanical  durability.  The  influence  of  processes  such  as  shot  peening,  induction  hardening,  and 
carburizing  on  fatigue  performance  depends  upon  the  magnitude,  distribution,  and  stability  of  re¬ 
sidual  stresses  developed  during  these  surface  processes.  Residual  stresses  may  also  develop  in 
service  if  localized  plastic  straining  occurs.  The  regions  of  compressive  residual  stress  are  well 
known  to  retard  fatigue  crack  growth  (or  arrest  cracks),  while  tensile  residual  stress  regions  pro¬ 
duce  the  opposite  effect.  However  the  overall  effect  of  residual  stresses  on  fatigue  life  is  much 
more  complex.  For  this  reason,  several  experimental  and  analytical  studies  have  been  conducted 
to  understand  the  behavior  of  residual  stresses  on  fatigue  life.  In  the  following  section,  a  brief  re¬ 
view  of  the  current  research  is  presented. 

4.2  Residual  Stresses:  Literature  Review 

The  influence  of  residual  stress  on  fatigue  crack  growth  has  been  approached  in  different 

ways.  Two  of  the  popular  ^’^p;  caches  are:  (a)  superposition  of  the  respective  stress  intensity  fac¬ 
tors  for  the  residual  stress  field  and  for  the  applied  stress,  and  (b)  a  crack  closure  model.  Nelson 
[1982]  presented  a  review  of  these  two  approaches  with  their  relative  advantages  and  drawbacks 
for  use  in  design  and  analysis.  In  this  paper,  only  the  Mode  I  crack  model  was  covered. 

In  the  superposition  of  stress  intensity  factors,  Kj  for  the  residual  stress  field  is  obtained  by 
loading  the  crack  faces  with  the  residual  stresses  which  exist  normal  to  the  plane  of  potential 
crack  growth  in  the  uncracked  body.  Such  residual  stress  intensity  factors,  have  been  com¬ 
puted  for  a  number  of  different  crack  types  and  crack  face  stress  distributions  (Underwood  and 
Throop,  1979;  Grant,  1978).  A  model  based  on  weight  functions  was  proposed  by  Parker  [1982]. 
By  using  weighting  functions  he  tried  to  predict  the  crack  surface  overlapping  effects  which  give 
rise  to  nonlinear  contact  conditions.  Further,  the  weighting  functions  were  used  in  an  iterative  pro¬ 
cess  to  correct  for  overlapping  effects.  The  weighting  function  theory  for  calculating  stress  inten¬ 
sity  factors  in  a  linear  elastic  cracked  body  was  introduced  by  Bueckner  [1970]. 

The  effective  stress  intensity  factors  are  then  taken  as  ^res  ** 

based  on  elastic  analysis.  This  is  the  major  drawback  of  this  type  of  analysis  as  it  lacks  the  ability 
to  account  for  the  influence  of  possible  changes  in  the  residual  stress  field  induced  by  service 
loadings  before  the  crack  starts  and  as  it  grows.  Further,  most  models  are  only  designed  to  predict 
growth  rate  under  constant  amplitude  loading. 

The  crack  closure  model  has  been  used  to  predict  load  sequence  and  mean  stress  effects  in 
specimens  with  crack-generated  residual  stress  fields.  The  significance  of  crack-generated  residu¬ 
al  stresses  was  pointed  out  by  Fiber  [1971],  who  observed  that  the  residual  tensile  displacements 
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left  in  the  wake  of  a  crack  would  cause  it  to  close  before  the  tensile  load  was  removed  and,  thus, 
generate  compressive  residual  stresses  in  the  wake.  The  crack  would  not  open  until  a  certain  level 
of  tensile  loading  was  applied  to  overcome  the  crack  opening  stress  resulting  from  the  compres¬ 
sive  residual  stresses.  By  correlating  crack  growth  rate  with  applied  mean  stress,  Elber  was  able  to 
provide  a  plausible  physical  explanation  for  retardation  caused  by  tensile  overloads.  This  method 
was  further  developed  by  other  researchers  (Dill  and  Saff  [1976],  Budiansky  and  Hutchinson 
[1978])  for  variable  amplitude  loadings.  The  primary  disadvantage  of  this  model  is  that  it  usually 
requires  an  elastic-plastic  finite  element  analysis,  which  is  computationally  cumbersome. 

A  recent  paper  by  Landgraf  et  al.  [1988]  formulated  a  stress-based  life  relationship  which 
incorporated  mean  stress  and  residual  stress  effects  to  identify  the  probable  failure  location  (sur¬ 
face  or  subsurface)  and  expected  fatigue  life.  In  this  model,  the  cyclic  properties  and  behavioral 
trends  developed  using  axial  specimens  of  steel  simulating  the  various  microstructure  found  in 
surface  treated  members  have  been  used  to  develop  a  comprehensive  procedure  for  analyzing  re¬ 
sidual  stress  effects  on  fatigue. 

Researches  have  also  been  conducted  on  the  optimization  of  surface  processing  methods 
and,  hence,  a  control  on  residual  stresses  which  would  provide  a  better  fatigue  life.  Fuchs  [1988] 
developed  a  method  of  analysis  for  finding  optimum  shot  peening  specifications  for  a  given  spec¬ 
imen  and  loading  conditions.  A  simplified  analysis  showed  that  the  optimum  shot  peening  intensi¬ 
ty  for  a  fillet  is  more  than  four  times  the  conventionally  specified  intensity.  The  effect  of  elevated 
temperature  on  shot  peened  steel  was  studied  by  Childs  [1988].  From  his  experiments,  he  con¬ 
cluded  that  a  significant  portion  of  the  residual  compressive  stress  is  lost  at  very  low  stresses.  He 
used  a  hole  drilling  method  to  measure  the  residual  stresses.  The  effect  of  residual  stresses  on  the 
morphology  of  fatigue  cracks  was  investigated  by  Knowles  et  a/.  [1991]  in  a  powder-metallurgy 
8090  aluminium  alloy.  Residual  stress  measurement  indicated  that  the  compressive  stress  levels 
developed  in  these  materials  was  significantly  greater  than  in  aluminum  with  conventional  metal¬ 
lurgy. 
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5.0  Summary  of  Work  Performed  in  FY  92 

Following  is  a  summary,  by  task,  of  the  research  performed  in  FY92. 

•  Task  1.  Improved  Experimental  Design  Methodology. 

Studies  have  been  performed  to  develop  optimal  design  strategies  fcu'  determining  the  tire- 
wheel  interface  loads  (normal  tractions)  for  landing  gear  wheels.  Although  an  accurate  represen¬ 
tation  of  the  tire-wheel  interface  loads  is  known  to  be  necessary  for  durability  and  damage  toler¬ 
ance  analysis  of  aircraft  wheels,  these  loads  are  difficult  to  obtain  by  either  experimental  or 
analytical  means.  Preliminary  research  has  developed  a  promising  combined  analytical/experi¬ 
mental  approach  to  overcome  these  difficulties.  This  methodology  requires  the  accumulation  of 
wheel  strain  gage  data  as  input  to  the  analytical  model.  However,  experience  to  date  indicates  that 
the  placement  and  orientation  of  the  strain  gages  are  crucial  to  successful  performance  of  the  al¬ 
gorithm.  Expanding  upon  initial  research  efforts,  UND  has  developed  a  methodology  for  optimal 
positioning  of  the  strain  gages  and  for  determination  of  the  minimum  number  of  gages  for  the 
analysis.  In  addition  to  the  Block  30  model  developed  under  the  previous  contract,  UND  has  de¬ 
veloped  a  finite  element  model  of  the  Block  40  F-16  main  landing  gear  wheel  in  support  of  this 
optimization  effort.  Specifications  for  the  required  number  of  strain  gages  and  their  optimal  posi¬ 
tion  for  the  Block  30  and  Block  40  F-16  main  landing  gear  wheels  were  provided  to  the  director¬ 
ate  at  WL/FIYMA  to  demonstrate  the  proposed  experimental  design  strategy. 

•  Task  II.  Tire* Wheel  Interface  Shear  IVactions. 

The  analytical/experimental  algorithm  has  been  extended  to  include  determination  of  tire- 
wheel  interface  shearing  tractions,  and  the  ANTWIL  finite  element  code  was  modified  to  include 
this  effect.  Complex  wheel  loadings  which  occur  in  braking,  spin-up,  cornering,  etc.,  produce 
both  normal  and  shear  tractions  at  the  tire-wheel  interface.  The  effect  of  the  shear  tractions  is 
thought  to  be  important  to  durability  and  damage  tolerance  of  the  wheel. 

•  Task  III.  Crack  Nucleation  Modeling  and  Residual  Stress  Effects  in  Aircraft  Wheels. 

Preliminary  research  into  models  of  crack  nucleation  (initiation),  including  the  effects  of  re¬ 
sidual  stresses,  has  been  initiated.  The  combination  of  cyclic  stresses  due  to  wheel  rotation  and 
corrosive  environmental  influences  can  cause  surface  cracks  to  form  which  further  enhance  stress 
concentrations  leading  to  crack  growth  and  ultimate  failure  of  the  wheel.  While  current  durability 
and  damage  tolerance  approaches  neglect  the  initiation  life  of  fatigue  cracks,  this  assumption  ap¬ 
pears  to  be  too  conservative  for  wheel  life  predictions  and  needs  to  be  reevaluated  in  the  context 
of  overall  wheel  life  prediction  methodology.  UND  will  report  on  the  most  promising  research  di¬ 
rections. 
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6.0  Planned  Tasks:  FY  93 


Given  below  are  the  tasks  planned  for  FY  93.  A  brief  description  including  motivating  argu¬ 
ments  is  given  with  each  task. 

•  Task  I.  Development  of  a  Three-Dimensional  Finite  Element  Model  of  the  Block  30  and 
Block  50  F-16  Main  Landing  Gear  Wheel 

A  fully  detailed  three-dimensional  solid  model  of  the  Block  30  and  Block  50  F-16  main 
landing  gear  wheel  will  be  developed  using  ARIES  [1992].  A  preliminary  solid  model  of  the 
Block  30  wheel  is  shown  in  Figure  27.  The  model  can  then  be  translated  into  PATRAN  [1988] 
format.  Using  this  solid  model,  finite  element  meshes  can  be  created  with  arbitrary  refinement. 
That  is,  the  mesh  can  be  modified  in  areas  of  high  stress  gradient,  such  as  the  bolt  hole  region,  to 
an  arbitrary  degree  of  refinement.  These  meshes  can  be  used  for  a  wide  variety  of  verification  and 
analysis  studies.  The  ability  to  assess  the  effects  of  bolt  torque  will  be  included  in  the  model. 
UND  will  provide  the  Directorate  with  the  full  three-dimensional  ARIES  solid  model,  as  well  as 
the  associated  finite  element  model. 


nCURE  27.  Model  of  F-16  Block  30  Main  Landing  Gear 

Wheel 
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*  Task  II.  Development  of  Methods  for  Applying  Tire* Wheel  Pressure  Distributions  to  the 
Three*Dimensional  Finite  Element  Model. 

The  tire-wheel  interface  pressure  distribution  determined  by  the  combined  analytical/exper¬ 
imental  procedure  developed  by  UND  is  essentially  an  analytical  representation  of  the  pressure. 
That  is,  the  Fortran  code,  ANTWIL,  yields  a  Fourier  series  representation  of  the  pressure  distribu¬ 
tion  in  the  coordinates  s  and  6  where  ^  is  a  coordinate  measuring  distance  along  the  wheel  flange 
and  6  is  the  angle  giving  circumferential  location.  Likewise,  any  pressure  distributions  obtained 
by  a  direct  measurement  technique,  recorded  digitally  or  continuously,  provide  essentially  a  con¬ 
tinuous  representation  of  the  interface  pressure  distribution.  In  either  case,  to  manually  input  the 
values  of  the  pressure  distribution  on  each  element  surface,  as  required  in  a  standard  three-dimen¬ 
sional  finite  element  code,  would  be  a  monumental  task  and  very  labor-intensive.  This  task  would 
have  to  be  repeated  for  every  simulation  with  a  different  pressure  distribution.  Since  many  com¬ 
puter  analyses  are  anticipated  with  different  pressure  distributions,  (from  different  tires,  different 
roll  tests,  etc.)  software  needs  to  be  developed  to  make  the  task  of  inputting  the  tire-wheel  pres¬ 
sure  distribution  into  a  three-dimetisional  finite  element  code  an  automated  process.  This  code 
should  be  able  to  take  the  required  coordinate,  nodal  point,  and  element  data  from  ARIES  and 
generate  nodal  work-equivalent  loads.  These  nodal  loads  should  be  output  in  a  fi'e  accessible  by 
the  user’s  finite  element  code  or  a  designated  commercially  available  code  such  as  NASTRAN. 
Upon  completion  of  this  task,  the  user,  given  the  ARIES  generated  three-dimensional  finite  ele¬ 
ment  mesh,  can  easily  perform  a  sequence  of  accurate  stress  analyses  for  any  measured  or  predict¬ 
ed  tire-wheel  interface  pressure  distribution.  UND  will  provide  the  Directorate  with  the  load 
generation  code  and  complete  documentation  including  sample  analyses.  The  load  generation 
code  will  interface  with  the  chosen  finite  element  code  given  the  finite  element  mesh  and  analyti¬ 
cally  or  experimentally  determined  interface  pressures. 

•  Task  III.  Verifleation  of  Experimentally  Determined  Pressure  Distributions. 

Use  of  a  pressure  transducer  offers  a  promising  technique  for  directly  measuring  tire-wheel 
interface  pressure  distributions.  Before  actual  application  to  the  complex  tire-wheel  problem,  the 
accuracy  of  any  experimental  determination  should  be  evaluated  by  controlled  experiments.  Thus, 
tests  need  to  be  designed  to  evaluate  the  capability  of  the  experimental  setup  to  accurately 
measure  interface  pressure  distributions.  The  experimental  design  should  be  such  that  the  pressure 
distribution  for  the  test  configuration  is  precisely  known.  In  addition,  the  tests  should  assess  the 
effect  on  the  accuracy  of  the  transducer  of  (a)  in-plane  stretching  of  the  transducer  that  will  likely 
occur  during  tire  mounting  and  (b)  bending  of  the  transducer  that  occurs  when  the  transducer  is 
mounted  in  the  bead  seat  region.  UND  will  provide  the  Landing  Gear  Section  of  the  Flight  Dy¬ 
namics  Directorate  with  the  design  of  appropriate  tests  and  analytical  solutions  with  which  to 
compare  the  experimental  results. 
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•  Task  IV.  Evaluation  of  Experimental  Setup  for  Tire- Wheel  Interface  Load  Determina¬ 
tion. 

Given  a  positive  evaluation  of  the  experimental  setup,  UND  will  input  the  tire-wheel  inter¬ 
face  pressure  distribution  obtained  from  the  pressure  transducers  during  wheel  roll  tests  into  the 
three-dimensional  Unite  element  model  and  calculate  the  resulting  strains  for  comparison  to  their 
experimental  counterpart.  This  will  provide  an  assessment  of  the  ability  of  the  experimental  set¬ 
up  to  perform  in  realistic  situations.  Upon  successful  verification,  UND  will  provide  the  Landing 
Gear  Section  of  the  Flight  Dynamics  Directorate  with  the  associated  evaluation  documentation. 

•  Task  V.  Verification  of  the  Proposed  Analytical/Experimental  Approach. 

Studies  will  be  performed  to  verify  the  proposed  analyticai/experimental  approach  for  tire- 
wheel  interface  load  determination  for  the  Block  30  F-16  main  landing  gear  wheels.  The  Landing 
Gear  Section  of  the  Flight  Dynamics  Directorate  at  WL/FIVMA  will  provide  UND  with  experi¬ 
mental  strain  survey  data  for  these  wheels  based  upon  the  experimental  design  specifications  pro¬ 
vided  in  Task  I.  Tire-wheel  interface  loads  will  be  extracted  from  the  experimentally  obtained 
strain  data  using  the  proposed  analytical/experimental  algorithm.  These  loads  will  then  be  reap¬ 
plied  to  the  three-dimensional  finite  element  model  of  the  wheel,  and  the  resulting  analytical 
strains  will  be  compared  to  their  experimental  counterpart.  In  addition,  these  load  profiles  will  be 
compared  with  those  obtained  from  the  laboratory  experiments.  Upon  successful  completion  of 
this  task,  UND  will  provide  the  Landing  Gear  Section  of  the  Flight  Dynamics  Directorate  with  the 
verification  documentation  and  the  derived  load  profiles. 

•  Task  VI.  Comparison  of  Tire- Wheel  Interface  Pressure  Distributions  for  the  Bias  Ply  and 
Radial  Tires  and  Their  Effect  on  the  Distribution  of  Wheel  Stresses. 

Bias  ply  and  radial  ply  tires  are  known  to  produce  different  stress  distributions  on  the  same 
wheels.  This  results  in  different  predictions  of  wheel  service  life.  WL/FTVMA  will  perform  roll 
tests  with  both  bias  ply  and  radial  ply  tires  and  the  strain  gage  results  will  be  made  available  to 
UND.  Using  the  combined  analytical/experimental  load  recovery  algorithm  implemented  in  the 
finite  element  code  ANTWIL,  UND  will  determine  the  interface  pressure  distributions.  These 
pressure  distributions  will  be  used  as  input  to  the  three-dimensional  finite  element  model  of  the 
wheel  developed  in  Task  I.  The  resulting  stress  distributions  in  the  wheel  will  be  determined  and 
compared.  The  comparison  will  include  assessment  of  the  effect  of  inflation  pressure,  vertical 
loads,  and  yaw.  UND  will  deliver  to  the  Landing  Gear  Section  of  the  Flight  E>ynamics  Directorate 
the  results  of  the  simulations  including  a  comparative  analysis. 

•  Task  VII.  Effect  of  Bolt  Torque  on  the  Distribution  of  Wheel  Stresses. 

The  magnitude  of  the  bolt  torque  is  known  to  have  a  significant  effect  on  the  wheel  service 
life.  This  undoubtedly  affects  the  stress  distribution  in  critical  areas.  Given  the  three-dimensional 
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finite  element  model  of  the  wheel  generated  in  Task  1  and  any  associated  software  and  the  ability 
to  easily  perform  repeat  analyses  (Task  V),  the  effects  of  the  bolt  torque  on  wheel  stress  distribu¬ 
tion  can  readily  be  studied.  UND  will  provide  the  Landing  Gear  Section  of  the  Flight  Dynamics 
Directorate  with  an  analysis  of  the  results  of  three-dimensional  finite  element  studies  performed 
to  assess  the  effects  of  different  bolt  torques. 

•  Task  VIII.  Residual  Stress  Measurements. 

UND  will  contract  with  Lambda  Research  of  Cincinnati,  Ohio  to  make  residual  stress  mea¬ 
surements  on  the  Block  30  F-16  main  landing  gear  wheel  at  several  critical  areas  and  at  several 
depths  for  inclusion  in  wheel  life  assessment.  A  specific  work  statement  will  be  written  for  Lamb¬ 
da  Research.  UND  will  keep  the  Landing  Gear  Section  of  the  Flight  Dynamics  Directorate  at  WU 
FIVMA  informed  of  the  progress  on  this  task. 

•  Task  IX.  Wheel  Fatigue  Life  Analysis. 

The  analytical/experimental  approach  for  determining  tire-wheel  interface  loads  will  be 
combined  with  fatigue  and  fracture  reliability  strategies  to  perform  wheel  life  analyses  of  aircraft 
wheels.  Both  stress-life  and  fracture  mechanics  approaches  (short  and  long  cracks)  will  be  con¬ 
sidered.  The  effect  of  residual  stresses  will  be  assessed.  WL/FIVMA  will  provide  roll-to-failure 
data  to  be  used  in  determining  possible  causes  of  wheel  failure.  WL/FIVMA  will  be  provided  ex¬ 
ample  analyses  for  prediction  of  wheel  life.  The  wheel  life  predictions  for  the  F-16  Block  30 
wheel  will  include  effects  of  stresses  due  to  bolt  torque,  inflation  pressure,  vertical  loads,  yaw,  and 
residual  stresses.  UND  will  also  investigate  possible  design  modifications  to  improve  wheel  life. 
A  subcontract  will  be  issued  to  obtain  required  residual  stress  measurements  for  the  F-16  Block 
30  main  landing  gear  wheel. 
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7.0  Closure 


This  report  has  deUneated  the  accomplishments  on  the  project  “Aircraft  Wheel  Life  Assess¬ 
ment”  during  FY92  and  has  reviewed  the  tasks  proposed  for  FY93.  Significant  progress  was  made 
in  developing  the  algorithm  for  back-calculating  tire-wheel  interface  pressures  with  the  goal  of 
making  this  a  reliable  component  of  the  overall  life  prediction  methodology.  The  ability  to  take 
experimental  strains,  back-calculate  the  tire-wheel  interface  pressure  and  apply  the  resulting  pres¬ 
sure  to  an  accurate,  three-dimensional  finite  element  model  of  the  wheel  will  be  totally  available 
sometime  during  FY93.  Complete  wheel  life  assessments,  including  a  comparison  of  bias  ply  and 
radial  ply  tires,  will  subsequently  be  performed. 

Additionally,  basic  literature  survey  work  into  the  difficult  questions  of  the  effects  of  resid¬ 
ual  stresses  on  the  wheel  life  assessment  problem  and  the  estimation  of  the  time  to  crack  initiation 
was  begun  in  FY92.  Identification  of  applicable  models  incorporating  these  effects  will  take  place 
in  FY93.  This  will  permit  the  life  estimation  methodology  to  be  considerably  improved  by  includ¬ 
ing  these  major  factors. 
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Appendix  A-  Finite  Element  Formulation 

This  appendix  will  present  the  concepts  of  finite  element  modeling  which  are  used  for  anal¬ 
ysis  of  axially  symmetric  bodies.  A  mathenuUical  model  for  axially  symmetric  bodies  is  de¬ 
scribed.  The  loading,  which  is  not  axially  symmetric,  is  also  discussed.  Loads  which  are 
asymmetric  are  solved  by  superposition  of  results  obtained  from  the  components  of  Fourier  se¬ 
ries. 

Introduction 

An  axisymmetric  solid  can  be  thought  of  as  a  solid  generated  by  revolving  an  area  in  the  n 
plane  around  the  z  axis,  as  indicated  in  Figure  A-1.  The  solid  body  which  is  generated  in  this 
manner  will  be  geometrically  axially  symmetric,  and  if  the  material  properties,  loads  and  support 
conditions  are  also  axially  symmetric  (i.e.,  the  variation  is  independent  of  the  d  direction),  then 
the  static  problem  will  have  displacements  and  stresses  which  are  independent  of  0.  This  essen¬ 
tially  means  that  the  problem  is  mathematically  two-dimensional,  where  displacements  are  de¬ 
fined  in  only  the  radial  and  axial  directions.  The  analysis  procedure  for  such  a  problem  is  similar 
to  that  of  plane  strain  or  plane  stress. 

In  most  cases  of  interest,  however,  the  solid  of  revolution  is  subjected  to  a  loading  which  is 
not  axially  symmetric,  or  the  material  property  may  not  be  isotropic.  In  either  case,  the  displace¬ 
ments  and  the  stresses  are  no  longer  two-dimensional  but  three-dimensional.  An  asymmetric  load¬ 
ing  can  be  represented  by  a  Fourier  series,  and  the  analysis  is  done  for  each  load  component.  Then 
using  the  principle  of  superposition,  the  original  problem  is  solved  by  adding  the  solutions  of  the 


r.  u 

FIGURE  A-1.  Axisymmetric  Body  with  IMangular  Elements. 
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component  problems.Thus,  a  three-dimensional  problem  is  broken  into  a  series  of  two-dimen¬ 
sional  problems.  In  most  cases  solving  the  problem  by  a  series  of  two-dimensional  components 
has  proved  to  be  fruitful  as  compared  to  a  full  three-dimensional  solution,  especially  in  terms  of 
memory  allocation  and  computational  time.  Besides,  three-dimensional  problems  are  difficult  to 
set  up  and  run. 

Formulation 

The  finite  element  method  for  linear  elastic  stress  analysis  can  be  derived  by  many  tech¬ 
niques.  The  approach  taken,  herein,  is  to  use  the  principal  of  virtual  displacements. 

The  strain-displacement  relationship  for  a  linear  elastic  material  can  be  written  as: 

€  =  LU,  (25) 

where  c  =  vector  of  strains,  L  =  matrix  operator  which  defines  the  strains  in  terms  of  the  displace¬ 
ments,  and  U  =  vector  of  displacements.  The  stress-strain  relation  is  given  by  the  following 
expression: 


<r  =  D(e-€j +a^,  (26) 

where  =  vector  of  initial  strains,  =  vector  of  residual  or  initial  stress,  and  D  =  elasticity 
matrix  which  depend  upon  the  properties  of  the  material. 

A  virtual  displacement  6U,  which  belongs  to  a  class  of  admissible  functions,  is  imposed  on 
the  system  and  the  external  work  is  equated  to  the  internal  work  and  integrated  over  the  entire  vol¬ 
ume  of  the  system.  This  relationship  is  expressed  as: 

J  ( )  d V  =  J  ( SV^h)  dV+  J  ( S\j\)  dA ,  (27) 

V  V  A 

where  b  =  vector  of  body  forces  per  unit  volume,  and  q  =  vector  of  surface  tractions  acting  over  a 
surface  A . 

The  formulation  of  the  basic  equations  can  be  done  at  the  element  level  where  the  equations 
are  summed  for  each  element  over  the  entire  continuum.  To  reduce  the  problem  to  finite  dimen¬ 
sions,  the  displacements  within  each  element  are  interpolated  from  nodal  displacements.  This  re¬ 
lation  is  given  by  the  following  expression: 


U  =  Na,  (28) 

where  a  =  vector  of  nodal  displacements,  and  N  =  matrix  of  interpolation  functions  or  shape 
functions.  Substituting  Equation  (28)  into  Equation  (25),  the  approximate  strain-displacement 
relationship  is  given  by: 


e  =  Ba, 


(29) 


where: 


B  =  LN,  (30) 

and,  £  -  approximate  strain  vector.  On  substituting  Equation  (28)  and  Equation  (30)  into  Equa¬ 
tion  (27),  the  virtual  work  expression  becomes: 

J(B5a)’'{D(Ba-e^) +<TjdV  =  J(N«a)’'bdV+f  (N5a)^q<M+5a^P  (31) 

V  V  A 

Rearranging  Equation  (31)  and  assuming  that  virtual  displacement,  5a,  is  not  equal  to  zero,  the 
above  expression  can  be  written  as: 


Ka  =  f ,  (32) 

where  K  -  stiffness  matrix  given  by, 

K  =  J(B'*'dB  )dV,  (33) 

V 

and  f  =  vector  of  loading  components  due  to  tte  body  forces,  traction  forces,  concentrated  forces 
and  residual  stresses  and  strains.  Further  it  can  be  expressed  as: 

f  =  J  (N^b)  dV+ J  (N^'q)  dA+j  (B^DeJ  dV+J  (BV„)  dV+P,  (34) 

V  X  V  V 

where  P  =  point  load  vector  acting  at  the  nodes. 

Axisymmetric  Body  with  Asymmetric  Loading 

As  alluded  to  earlier,  axisymmetric  solids  may  have  loads  which  may  or  may  not  be  axially 
symmetric.  The  situation  in  which  loading  is  asymmetric  can  be  solved  by  decomposing  the  load¬ 
ing  function  using  a  Fourier  series  in  the  ^-direction  and  superposing  the  results.  The  body  will 
have  deformations  in  the  r,  $  and  z  directions.  In  fact,  the  symmetric  case  can  be  seen  as  a  special 
case  of  the  above,  by  having  only  the  zeroth  component  of  the  Fourier  series  (Figure  A-2a)  where 
the  displacements  are  only  in  the  r  and  z  direction.  Thus,  a  mathematical  model  for  the  most  gen¬ 
eral  case  can  be  developed  which  incorporates  the  special  case  of  symmetric  loading. 

Let  the  loading  on  the  body  be  expressed  as  a  Fourier  series,  which  can  be  written  as  [Cook, 
1977]: 


45 


(35) 


[Qr  .Qj  .Sj  .T  ]  =  ^V„cCos  (nO)  ind) 

n  n 

[Q^S^J  =  '^L\s5inind)-*-'^V‘„cCosind)  (36) 

n  n 

where  Q’s  and  S’s  are  vectors  of  body  forces  and  surface  tractions  in  the  r,  0  and  z  directions 
and  T  =  vector  of  temperatures.  L  contains  the  amplitudes  of  the  vectors  of  body  forces  and 
surface  tractions.  and  are  the  cosine  and  sine  amplitudes  and  represent  a  state  of 
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symmetry  with  respect  to  the  plane  0=0,  where  as  and  are  the  sine  and  cosine  ampli¬ 
tudes  for  the  antisymmetric  case.  These  equations  can  be  expressed  as: 


(37) 

K  =  (Ci,  •  ■'‘•.1  <3*) 

1.;,=  le;,,  .St,.  (39) 

(«) 

It  can  be  shown  [Cook,  1977]  that  under  the  above  loadings  the  displacements  have  the  fol¬ 
lowing  form: 

u  =  ^u^cos  (n^)-i-2^u“sin  (n0)  (41) 

n  n 

V  =  ^v^sin  (n^) -^v®cos  (n^)  (42) 

n  n 

w  =  ^w^cos  (n0)+2^w“sin  (nd)  (43) 

n  n 


where  u  =  radial  displacement,  v  s  circumferential  displacement  and  w  =  axial  displacement. 

The  strain-displacement  or  the  kinematic  relation  which  is  given  in  Equation  (26)  is  ex¬ 
pressed  here  in  cylindrical  coordinates: 

'd/idr)  0  0 

l/r  d/(rdd)  0 

0  0  9/  (dz) 

9/(9z)  0  d/(9r) 

d/(rdff)  (d/iBr)-l/r)  0 

0  d/ (dz)  d/(rd0) 

The  stress-strain  relationship  is  given  by  Equation  (26)  ,  o-  =  D(e-  e^)  -i-  (r^  where  D,  the 
constitutive  matrix,  is  expressed  for  isotropic  material  as  [Timoshenko,  1970]: 
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(45) 


1  /  /  0  0  0 

/  1  /  0  0  0 

/  /  1  0  0  0 

0  0  0  g  0  0 

0  0  0  0  G  0 

_0  0  0  0  0  GJ 

where  v  =  poisson  ratio,  E  =  Young’s  modulus  of  elasticity  for  the  material  and  G  =  shear 
modulus  of  elasticity.  Also  /  and  g  can  be  expressed  in  terms  of  Poisson’s  ratio  and  Young’s 
modulus  as: 


D 


( 1  -  i.)  £ 


(1  +  t^)  (l-2i/) 


/  = 


V 


^  {\-2v) 
2(l-p) 


(46) 


The-  contribution  of  the  harmonic  to  the  strain-displacement  relationship  is  obtained  by 
substituting  Equation  (41)  through  Equation  (43)  into  Equation  (44).  This  relationship  can  be  ex¬ 
pressed  as  €„  =  B„a„,  or  more  specifically  for  the  symmetric  component  by: 


N,  cos  (nff) 

*  t  Z 

0 

0  ..' 

r  ■ 

^1 

nN. 

“  *• 

—  cos  (nd) 

— -cos  (nff) 

0 

“i 

r 

0 

r 

0 

N.  ,cos{nd)  ... 

: 

yrz 

N^  jCos  (nd) 
-nN, 

0 

(yv,  ^-^^)sin(/i^) 

IV,  ^cos(n0)  ... 

yre 

— ^sin  (nff) 

0 

... 

w 

0 

N,  ,sin(n0) 

—nN, 

— - — sin(n^)  ... 

1.2,3 . 8 

(47) 


It  should  be  noted  that  for  antisymmetric  harmonics,  the  B  matrix  is  obtained  by  inter¬ 
changing  the  sine  and  cosine  terms,  and  the  last  two  rows  are  multiplied  by  negative  sign.  The 
stiffness  matrix  is  obtained  by  substituting  Equation  (47)  in  Equation  (34).  One  also  finds  that  the 
stiffness  matrix  obtained  for  the  symmetric  and  antisymmetric  cases  are  identical.  This  is  due  to 
the  nature  of  displacements  given  by  Equations  (41),  (42)  and  (43).  The  stiffness  matrix  obtained 
above  is  only  a  function  of  n.  Thus,  the  assembly  of  all  individual  stiffness  matrices  results  in  a 
system  of  uncoupled  equations  which  can  be  solved  independently  for  each  n.  The  solutions  ob¬ 
tained  from  the  analyses  are  then  superposed  to  get  the  overall  result. 
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Appendix  B  -  Surface  TVactions  and  Load  Distributions 

This  chapter  discusses  the  theory  required  for  the  particular  implementation  of  distributed 
surface  tractions  employed  in  this  research.  Strain  and  stress  computations  at  any  point  on  the  sur¬ 
face  are  also  discussed. 

Discretization  of  Surface  TVactions 

In  general,  the  load  vector  f  given  in  Eq.  (34)  has  contributions  from  body  forces,  surface 
tractions,  residual  stresses  and  strains,  and  concentrated  forces.  One  can  use  any  shape  functions 
to  interpolate  the  distributed  loads,  but  in  most  cases  the  shape  functions  used  are  the  same  as 
those  used  to  develop  the  stiffness  matrix  and  the  loads  obtained  are  called  consistent  nodal  loads. 
In  this  section,  only  the  surface  tractions  for  the  eight  noded  element  are  discussed.  Quadratic 
shape  functions  are  used  to  interpolate  the  load  distribution  and  hence  the  loads  obtained  are  con¬ 
sistent  [Bickford,  1990  and  Cook,  1989]. 

Let  q  (j)  represent  the  surface  traction  vector  on  side  three  of  the  element  shown  in  Figure 
B-1.  The  vector  qf^)  has  components  and  qyis) .  The  variable  s  represents  distance 

along  the  surface  defined  by  nodes  3, 4,  and  7.  The  contribution  to  the  load  vector  from  the  load¬ 
ing  on  this  element  is  then  given  as: 


f^is)  =  jN^q(5)</5 

5 


(48) 


FIGURE  B-1.  Eight-Noded  Isoparametric  Element  with  a  Load  of  q(s). 
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where  is  a  ,6  X  1  vector  of  consistent  nodal  loads  and  N  is  a  2  x  6 matrix  of  element  shape 
functions.  To  employ  gauss  integration,  a  natural  coordinate  system,  17,  is  introduced.  In  the 
natural  coordinate  system,  the  element  is  a  2  x  2  square  where  ^  and  rj  each  run  from  -1  to  +1. 
Using  the  element  shape  functions  to  also  define  the  mapping  between  physical  coordinates  and 
natural  coordinates,  x  and  y  can  be  expressed  as  [Bickford,  1990]: 

i  ' 

and: 


i  ' 

The  index  i  is  summed  over  nodes  3, 4  and  7,  the  nodes  on  side  3.  Note,  also,  that  side  3  is 
defined  in  the  natural  coordinate  system  by  77  =  1 .  Further  the  derivatives  of  x  and  y  with  re¬ 
spect  to  ^  can  be  given  by  the  following  expressions: 


and: 


^  _  a 
di 


=  A^+B 


(51) 

(52) 


=  C^+D 

where  A  =  {x-^-lx-j^ x^) ,  B  =  {x^-x^)/2,  C  =  ( J3 - 2yy ^4)  and  D  =  (>'3-y4)/2. 
In  the  above,  (x^ ,  y^)  are  the  coordinates  of  node  i  (either  3, 4  or  7).  Sinularly,  the  differential 
length  ds  can  be  given  as: 


(53) 

(54) 


ds  =  Jdi 


(55) 


where: 


J  = 


2  2 

(-) 

vde  ^d^  J 


(56) 
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Furthermore  the  differential  length  ds  can  be  integrated  to  get  5,  the  distance  along  the  surface 
from  node  3,  that  is, 

S 

s  =  jds 
0 

( 

^jjdi  (57) 

-1 

Alsoi,  J  can  be  expressed  in  terms  of  A,  B,  C  and  D  as: 

I  /2 

J=  [(A^  +  C^)^  +  2(AB  +  CD)^+ (B^  +  D^)l  (58) 

It  can  be  seen  from  Figure  B-2  that: 

q^ds  =  pdy  (59) 

and 

q^ds  =  pdx  (60) 

Thus,  the  x  and  y  components  of  the  vector  f  can  be  written  as: 

/"ix  =  Jb  (^) 
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and: 


f/y  =  \pis)N,dx  (62) 

5 

where  again  /  is  3, 4,  or  7.  Using  Equations  (51)  through  (58),  these  integrals  can  be  evaluated  in 
the  natural  coordinate  system  as: 

I 

(Ci+D)di  (63) 

-1 

and: 

1 

/,;  =  (A^+B)di  (64) 

-I 

Computation  of  Strains  on  the  Surface  of  an  Element 

In  general,  strains  are  obtained  from  the  discretized  displacement  field  through  use  of  the 
strain-displacement  relation,  c  =  Bu,  where  c  =  vector  of  strainr  B  =  strain-displacement  ma¬ 
trix  and  u  =  vector  of  nodal  displacements  [Cook,  1989;  Timosiicnko  and  Goodier,  1970].  The 
important  question  which  arises  here  is  which  strains  should  be  computed  to  get  maximum  accu¬ 
racy.  It  can  be  shown  that  for  isoparametric  elements,  strains  are  more  accurate  at  the  Gauss 
points  of  the  quadrature  rule  one  order  less  than  that  required  for  full  integration  of  the  element 
stiffness  matrix  (i.e.,  the  so-called  Barlow  points).  However,  in  most  cases  the  points  of  interest 
are  not  the  Barlow  points.  In  such  cases,  the  strains  are  interpolated  or  extrapolated  to  any  point  of 
interest  from  the  Barlow  points  [Barlow,  1976]. 

Extrapolation/Interpolation  of  Strain  Tensor 

Consider  the  eight  noded  isoparametric  element  depicted  in  Figure  B-3.  Generally,  the  order 
of  the  Gauss  quadrature  rule  used  to  obtain  the  element  stiffness  matrix  is  3  x  3.  This  means  that 
the  strains  are  more  accurate  at  the  2  x  2  Gauss  points.  Furthermore,  it  can  be  deduced  from  the 
shape  functions  and  the  strain-displacement  relationship  that  displacements  are  quadratic  in  na¬ 
ture  and  strains  are  linearly  related.  Hence,  the  strains  at  any  point  on  the  surface  can  be  linearly 
interpolated  from  the  Barlow  points  A,  B,  C,  and  D. 

Consider  the  case  where  strains  are  to  be  derived  at  point  P  (see  Figure  B-2).  This  is  done  in 
two  stages.  At  the  first  stage,  strains  are  computed  at  the  four  Barlow  points  from  the  strain-dis¬ 
placement  relationship.  These  points  are  represented  in  the  r  -  s  coordinate  system.  The  r  -  5  co¬ 
ordinate  system  is  linearly  related  to  the  17  coordinate  system  by  the  following  expression: 
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r  =  75^  and  s  -  Srj  (65) 

In  second  stage,  the  strains  computed  at  the  Barlow  points  are  extrapolated  to  obtain  the 
strains,  Cp,  at  point  P  via: 


4 

«/.  =  X 

i=  1 

where  e,  =  vector  of  strain  at  the  Barlow  points  and  =  bilinear  shape  function  given  in  the  r-  s 
coordinate  system  for  the  four  Barlow  points  A,  B,  C,  and  D.  The  shape  functions  are  locally 
described  as. 


yv.  =  l(l±r)(l±r)  (67) 

where  r  and  s  are  the  local  coordinates  of  point  P. 

Computation  of  Strains  at  an  Arbitrary  Orientation  on  the  Surface 

Consider  the  situation  where  the  direct  strain  has  to  be  computed  in  the  direction  of  unit 
vector  g  at  an  angle  <p  with  respect  to  the  unit  base  vector  on  the  surface  of  an  element  (sec 
Figure  B-4).  At  any  point  on  the  surface  a  unit  normal  vector,  n,  and  a  unit  tangential  vector,  e,, 
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FIGURE  B>3.  Eight'Noded  Element  A,  B,  C,  and  D  are  Barlow  Points. 
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FIGURE  B-4.  (a)  Axtsymmetric  Body  (b)  Direction  of  g  at  Which  Strains  arc  Compatcd. 


are  computed  from  the  shape  functions  and  the  coordinates  of  nodes  y,  k  and  P  (Figure  B-5).  It 
has  to  be  noted  that  j  and  k  are  nodes  of  an  element  where  P  is  any  point  on  the  surface.  The  co¬ 
ordinates  of  P  are  computed  from  Equations  (49)  and  (50).  Using  Equations  (51)  through  (56), 
the  unit  normal  and  tangential  vectors  are  expressed  as  (see  Figure  B-4): 


and: 


en(^ 


1  fClx  dy  . 

di 


=  /(^)e^  +  m(^)ej 


*-<«  =  i 

di 


r  dy  dx  , 


=  -m(^)e,  +  /(Oe, 


(68) 

(69) 


I 

(70) 


(71) 
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FIGURE  B<5.  Side  of  an  Eight-Noded  Isoparametric  Element  Where  i)  =  1. 
Point  P  is  Where  Strains  arc  Compnted. 


where  and  are  unit  vectors  described  in  r  and  z  direction,  respectively.  Furthermore  I 
and  m  are  described  by  : 


no 


Ai+B 

J 


and: 


(72) 


m(^ 


C^+D 

J 


(73) 


where  A,  B,  C,  and  D  have  the  same  expression  as  defined  in  the  beginning  of  Appendix  B.  From 
Figure  B-4b  it  is  noted  that  the  vector  g  is  at  an  angle  <p  with  respect  to  in  the  -  e,  plane. 
Thus: 


g  =  cos  ( <p)  e^+  sin  ( <p)  e,  (74) 

Using  Equation  (71)  and  substituting  in  above  equation,  we  have: 

g  =  -m(^)sin(^)e^  +  cos(^)e^+/(^sin(^)ej  (75) 

or  in  matrix  notation: 
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(76) 


*  -  sin(<p)  cos(<p)  /(^sin(^)] 

where  again  is  a  unit  vector  in  the  6  direction.  Furthermore,  it  can  be  shown  that  the  required 
strain  vector  A  can  be  expressed  as: 


A  =  g^'eg  (77) 

where  c  is  the  strain  tensor  described  in  tt^  cylindrical  coordinate  system  by  the  following 
matrix. 


111  hi 

*'•2  2 

yre  yz0 
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(78) 
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